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(i) 


ABSTRACT 


A  recent  experiment  has  shown  that  a  chemical  system  under¬ 
going  a  reaction  may  exhibit  spatial  structure.  In  this  experiment  the 
reaction  mixture  gradually  arranged  itself  into  alternating  layers  of  red 
and  blue.  In  the  blue  layers  the  concentration  of  one  of  the  species 
present  was  higher  than  its  concentration  in  the  red  layers.  This  paper 
will  discuss  a  possible  theoretical  explanation  for  the  appearance  of 
structure  in  chemical  systems.  It  is  shown  that  this  may  happen  if  the 
homogeneous  steady-state  solution  to  the  chemical  rate  equations  (including 
diffusion)  can  become  unstable  with  respect  to  spontaneous  spatially  depend¬ 
ent  fluctuations. 

The  first  chapter  of  this  paper  is  an  introduction,  including  a 
historical  summary.  The  second  chapter  defines  the  concepts  of  stability, 
marginal  stability,  and  instability  for  homogeneous  steady-state  solutions. 
The  two  types  of  boundary  conditions  are  discussed  and  a  criterion  is 
developed  for  each  type  of  boundary  condition  to  determine  if  a  homogeneous 
steady-state  solution  is  stable  or  not.  The  criterion  is  valid  for  any 
dimension  and  shape  of  system,  provided  that  the  surface  of  the  system  is 
smooth. 


Next  a  specific  example  of  a  chemical  system  which  can  become 
unstable  with  respect  to  spatially  dependent  fluctuations  is  presented^. 
This  is  followed  by  an  example  which  illustrates  that  the  behavior  of  a 
system  can  differ  depending  on  which  type  of  boundary  condition  is  imposed 
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(ii) 
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on  it.  Following  this  a  simple  model  developed  by  Lefever  *  ’  is 

discussed  and  some  numerical  results  are  presented. 

The  fourth  chapter  is  concerned  with  the  concept  of  multiple- 
steady-state  solutions  to  the  chemical  rate  equations.  A  simple  mechanism 
is  presented  and  the  stability  of  all  the  homogeneous  steady-state  solutions 
is  studied  in  detail  using  Lefever's  model.  Then  it  is  shown  that  a  homogen¬ 
eous  steady-state  solution  can  be  stable  simultaneously  with  an  inhomogeneous 
steady-state  solution. 

The  conclusion  of  this  paper  will  discuss  in  more  detail  the  labor¬ 
atory  experiment  which  confirms  that  spatial  structure  can  exist  in  chemical 
systems . 
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CHAPTER  1 

INTRODUCTION 


It  is  common  to  think  of  chemical  reactions  occurring  in  a 

spatially  homogeneous  medium  -  that  is,  the  concentrations  of  the  various 

chemical  species  are  distributed  uniformly  throughout  the  reaction  vessel. 

At  thermodynamic  equilibrium  and  in  near-equilibrium  situations  this  must 
14 

be  the  case.  This  is  true  because  any  spontaneous  fluctuation  which 

takes  the  system  away  from  the  homogeneous  state  must  die  out  in  time. 

Thus  the  system  cannot  evolve  away  from  the  homogeneous  situation.  However, 
in  systems  far  from  equilibrium  this  may  not  be  the  case.  It  might  be  that 
fluctuations  could  grow  in  time  and  the  system  would  evolve  away  from  the 
homogeneous  state. 


Mathematically,  the  study  of  fluctuations  can  be  linked  to  the 

study  of  the  stability  of  differential  equations.  For  chemical  systems, 

the  appropriate  differential  equations  are  the  chemical  rate  equations 

including  the  effect  of  diffusion.  This  problem  was  first  studied  by 
21 

Turing  in  1952.  In  his  paper,  he  presented  several  reaction  mechanisms 

which  could  lead  to  a  structured  system.  However,  after  the  publication  of 

this  paper,  little  work  was  done  in  the  chemistry  field  for  a  long  time. 

But  the  general  area  of  stability  was  investigated  more  by  people  in  hydro- 
2  20 

dynamics  .  Then  in  1966  Scriven  and  Gmitro  published  a  paper  recalling 
Turing's  work  and  linking  chemical  systems  to  the  type  of  work  which  had 
been  done  in  hydrodynamics.  In  1967  and  1968  Prigogine^  *  ^  and  several 
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other  people  connected  with  him  *  published  some  articles  along  the 
same  lines  as  Scriven  and  Gmitro’s,  but  including  an  explanation  of  the 
phenomena  of  dissipative  instability  in  terms  of  non-equilibrium  thermo¬ 
dynamics.  Since  then  several  more  articles  and  a  book  have  been  published. 

This  paper  continues  the  study  of  spatial  structure  in  chemical 
systems.  First  a  mathematical  criterion  for  determining  stability  or  lack 
of  it  is  developed.  This  criterion  is  then  applied  to  a  mechanism  to  give 
a  straightforward  example  of  a  dissipative  instability.  Then  the  question 
of  the  boundary  conditions  imposed  on  the  system  is  discussed  at  length. 

This  problem  seems  to  have  been  ignored  by  most  authors.  However,  very 

recent  works  by  Ross  and  Ortoleva  »  have  also  discussed  boundary  behavior, 
although  their  treatment  is  very  different  from  the  approach  taken  in  this 
paper.  Following  this,  an  example  of  the  phenomena  of  multiple  steady- 
states  is  presented  and  the  stability  of  each  state  is  examined.  A  further 
extension  of  the  concept  of  multiple  steady-states  is  then  discussed.  This 
part  of  the  paper  was  motivated  by  an  unsolved  problem  presented  by  Edelstein^. 
This  extension  of  multiple  steady-states  shows  for  a  simplified  model  an 
answer  to  Edelstein's  problem.  The  general  problem,  as  yet,  remains  unsolved. 
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CHAPTER  2 

MATHEMATICAL  PRELIMINARIES 


For  the  rest  of  this  work  we  shall  be  concerned  about  chemical 
systems  in  which  the  reactions  are  assumed  to  be  completely  known.  For 
example,  the  system  could  be  governed  by  the  following  reaction  sequence 
or  mechanism: 


A 


X  — >  Y 

Y  +  B  — >  C  . 


The  overall  reaction  is  A  +  B  C  .  It  will  be  assumed  that  the  concen¬ 
trations  of  A  ,  B  ,  and  C  are  fixed  and  constant  in  time.  Then  the 
system  of  differential  equations  (chemical  rate  equations)  associated  with 
this  system  is: 


dX 

dt 


k  A  -  k2X 


dY 

dt 


k2X  -  k3BY 


Capital  letters  are  used  in  the  chemical  rate  equations  to  denote  the 
concentration  of  the  corresponding  substance.  This  system  has  one  time 
independent  or  "steady-state"  solution  obtained  by  setting 


J'.'.s  ;'] 


■ 


■ 


and  solving  the  resulting  algebraic  equations.  In  this  case  the  steady- 
state  solution  is: 


X 

o 


k  A 

Y  =  — — 
o  k„B 


In  general,  for  a  system  with  N  time  dependent  species 
(intermediates)  we  have: 

dX 

=  fn(X1,X2,-“,Xn)  ;  n  =  1 ,  •  •  •  ,N  (2.1) 

For  the  present  we  will  assume  there  exists  one  steady-state  solution, 

v  y  •  •  •  y 

10*  20*  *  NO  * 

Let  us  now  consider  solutions  that  have  been  perturbed  from 
the  steady-state  solution.  These  can  be  written  as: 

xn(t)  =  xn0  +  ’  n  =  (2.2) 

The  ljj  (t)  represent  fluctuations  away  from  the  steady-state  solution. 

These  could  represent  Brownian  motion  or  some  other  spontaneous  disturbance 

in  the  system.  The  ^  (0)  are  assumed  to  be  small  in  comparison  with 

X  _  .  Physically,  we  would  expect  that  a  small  disturbance  in  the  system 
nO 

would  not  drastically  alter  the  system.  That  is,  we  expect  that  the 

perturbed  solutions  will  remain  near  the  steady-state  solutions  as  time 

increases.  This  will  be  the  case  if  each  lb  (t)  +  0  as  t  ->  00  .  Then  we 

n 

say  that  the  steady-state  solution  X^q  ;  n  =!,•••, N  is  stable.  If  any 


•  v* 


I 

. 


5 


ip  (t)  becomes  large  in  comparison  to  ip  (0)  ,  we  will  say  that  the 


n 


n 


steady-state  solution  is  unstable.  If  none  of  the  ip  (t)  become  large 

n 

in  comparison  to  ip  (0)  and  at  least  one  of  the  ip  (t)  does  not  tend  to 

n  n 

zero,  then  we  will  say  that  the  steady-state  solution  is  marginally  stable. 


Let  us  now  substitute  equation  (2.2)  into  (2.1).  This  gives 


dip 


n 


dt 


fn^X10  +  ^1^  *  X20  +  ^2^  ’ 


’  XN0  + 


(2.3) 


n  =  1 ,  • • • ,N 


We  will  expand  the  right  hand  side  in  a  Taylor  series  and  retain  only  terms 

of  order  ip  (t)  .  We  obtain: 
n 


dip 


n 


N  /  3f 


^ 


n 


dt  nv  10’  20 


-  ,  3X . 

n=l  \  j 


X1  X10’”',XN  \o' 


dip  N  /  3  f 

—  =  y  — 

/J3X. 

1=1\  J 


dt 


X10’X20’ 


■^0 


(2.4) 


n  =  1 ,  •  •  •  ,N 


We  have  used  here  the  fact  that  the  X  _  are  the  steady-state  solution. 

nO 

Equations  (2.4)  are  now  a  system  of  linear  differential  equations  for  the 

fluctuations.  Provided  that  it  was  valid  to  retain  only  terms  of  order 

ip  (t)  (linearization  of  equations  (2.1)),  we  can  now  solve  (2.4)  and  find 
n 

the  fluctuation  explicitly  as  functions  of  time  and  from  this  determine 
the  stability  of  the  steady-state  solution.  The  validity  of  linearization 


' 
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of  equation  (2.1)  is  determined  by  the  form  of  the  functions  f  .  We 
will  make  use  of  the  following  theorem^. 

Given  a  system  of  differential  equations: 
dX 

— =  fn(X1>---,XN)  ;  n  =  1 ,  •  •  •  ,N  (2.5) 


such  that  each  f  can  be  written  in  this  form: 

n 


f 

n 


N 


I 

j-i 


a  . 
nJ 


X,  + 
1 


;n(Xl> 


n  =  1 ,  •  *  •  ,N 


(2.6a) 


where  the  a  are  constants,  and  such  that: 
n  . 

J 


V 


+ 


0  as 


+  XN*° 


n  =  1 


>  N 

(2.6b) 


where  t,  is  the  half  life  of  the  overall  reaction, 

then  if  the  corresponding  linearized  system  of  differential  equations  is 
stable,  system  (2.5)  is  stable  or  marginally  stable.  If  the  correspond¬ 
ing  linearized  system  of  differential  equations  is  unstable,  system  (2.5) 
is  unstable. 


We  will  find  that  in  all  the  cases  we  wish  to  consider  conditions 
(2.6a)  and  (2.6b)  will  be  satisfied. 

We  are  now  in  a  position  to  solve  equation  (2.4)  explicitly. 

Let  ip  be  the  column  vector  whose  nth  element  is  ip  and  let  A  be 
—  n  = 


the  matrix  defined  by 


■ 
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xrxio-->w 


(2.7) 


Then  (2. A)  can  be  rewritten  in  matrix  notation  as: 

dip 

—  =  .  (2.8) 

Once  it  is  in  this  form  it  is  easy  to  show  that  the  solution  is: 

At 

ip(t)  =  ip^e 


Further  each  fluctuation  ^  may  be  expressed  as: 


ip  =  )■  c  .e  J  ;  n  =  1 ,  •  •  •  ,N 

n  j=l  nJ 


(2.9) 


where  c  is  some  constant  amplitude  and  the  w.  are  the  eigenvalue 
nJ  3 

of  the  matrix  A  .  Clearly  if  the  real  part  of  each  w^  is  less  than 

zero,  each  fluctuation  ip  will  approach  zero  as  time  increases.  Then 

the  steady-state  solution  will  be  stable.  However,  if  just  one  of  the 

eigenvalues  has  a  positive  real  part,  the  perturbation  will  increase  in 

time  and  the  steady-state  solution  will  be  unstable.  If  one  of  the  w 

should  have  a  real  part  equal  to  zero,  this  would  correspond  to  the  case 

of  marginal  stability.  However,  examples  can  be  constructed  in  which  the 

linearized  system  of  equations  is  marginally  stable  and  the  corresponding 

nonlinear  system  is  not.  So  if  we  find  an  w.  =  0  we  cannot  conclude 

3 

that  the  steady-state  solution  to  a  nonlinear  set  of  chemical  rate  equations 
is  marginally  stable.  However,  a  positive  real  part  of  one  of  the  w^ 
would  still  indicate  instability.  Thus  we  see  that  determining  the  stability 


:  ’• 
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of  the  steady-state  solution  for  a  system  of  chemical  rate  equations  is 
equivalent  to  determining  if  the  eigenvalues  of  the  matrix  A  corres¬ 
ponding  to  the  system  have  any  positive  real  parts. 

Clearly  A  will  depend  on  the  rate  constants,  and  on  the  con¬ 
centrations  of  the  reactants  and  possibly  the  products.  So  the  eigenvalues 
of  A  will  be  dependent  on  the  same  quantities  and  it  may  turn  out  that 
the  system  becomes  unstable  for  some  values  of  the  rate  constants  and 
concentrations.  This  means  that  in  the  instability  region  if  a  fluctuation 
occurs,  the  solution  to  the  kinetic  equations  will  be: 


X 

n 


Xn0  +  V0 


n  =  1 ,  •  •  •  ,N 


and  this  solution  will  go  away  from  the  steady-state  solution.  That  is, 
the  concentration  of  the  will  become  time  dependent.  An  example  of 
this  will  be  considered  later. 


Until  now  we  have  not  considered  any  part  of  our  problem  to  be 
spatially  dependent.  Let  us  now  allow  diffusion  of  the  intermediate  species 
in  the  system.  Using  Fick’s  law  of  diffusion  the  chemical  rate  equations 
(2.1)  now  become: 


9X 

it-  fn(Xl-X2 


,XJ  +  D  V  X 
N  n  n 


(2.10) 


n  =  1 


,N 


D  is  the  diffusion  constant  of  the  nth  intermediate  species, 
n 


. 

' 


■ 


■ 
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It  should  be  noted  that  we  are  now  considering  a  system  of 
coupled  partial  differential  equations.  It  will  be  necessary  to  impose 
boundary  conditions  on  this  system  of  equations.  In  a  system  with 
diffusion  two  different  types  of  boundary  conditions  may  be  applicable 
depending  on  the  experimental  situation^.  One  is  to  have  fixed  values 
of  the  concentrations  at  the  boundaries.  An  example  of  this  is: 

Xn(£,t)  =  X^q  on  3S  ;  n  =  1,*‘*,N  (2.11) 

where  9S  is  the  boundary  of  the  region  S  in  which  the  system  under 
consideration  is  defined.  This  type  of  boundary  condition  is  applicable 
if  there  is  diffusion  through  the  boundary,  for  example,  with  a  membrane. 
This  kind  of  boundary  condition  will  be  referred  to  as  the  fixed  boundary 
condition.  However,  if  the  body  which  composes  the  boundary  is  insoluble 
in  the  fluid  in  the  system  and  so  diffusion  cannot  take  place  across  the 
boundary,  then  the  appropriate  boundary  condition  is: 

X 

— =0  at  the  boundary;  n  =  1,***,N  (2.12) 

where  v  is  the  unit  outward  normal  to  the  surface  of  S  .  We  would 
expect  this  type  of  boundary  conditions  to  be  applicable  in  most  experi¬ 
mental  situations.  We  will  refer  to  this  type  of  boundary  conditions  as 
derivative  boundary  conditions. 

We  will  assume  that  equation  (2.10)  has  a  "steady-state  homogen¬ 
eous"  solution  obtained  by  solving: 


\ 


10 


fn(xi,x2, ••• ,XN)  =0  ;  n  =  1 , • • • ,N 

We  will  designate  this  solution  by: 


(2.13) 


xnQ  ;  n  =  1 , *  *  * ,N  .  (2.14) 

We  again  wish  to  investigate  the  stability  of  this  solution. 

As  before  we  write  the  perturbed  solution: 

x  =  x  n  +  ip  (r , t)  ;  n  =  1,***,N  . 

n  nO  n  —  ’ 

Note  that  the  fluctuations  are  now  functions  of  _r  as  well  as 
Substituting  (2.15)  into  (2.11)  gives: 

'  (£>t)  _  _  2- 

~ 5^“  -  fn<X10+V->  W  +  V 

n  =  1 ,  •  *  •  ,N 


(2.15) 

t  . 


(2.16) 


Again  we  expand  in  a  Taylor  series  keeping  only  terms  linear  in  ip^(rft)  . 
We  obtain: 


dip  (r,t)  N 

dt 

n  ~  -  y 

n 

3t 

3X . 

j=i 

J 

X  =X  • • • j 

1  10* 

ip.  +  D  V2ip 
J  n  n 


(2.17) 


We  will  now  need  to  examine  what  effect  each  type  of  boundary  condition 

will  have  on  the  perturbations.  First  we  will  consider  the  fixed  boundary 

condition.  Since  the  x  _  ,  n  =  1***,N  are  solutions  we  must  have  that: 

nO 


\ 
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X  =  X  A  +  ip  (r,t)  =  X  A  on  9S  ;  n  =  1,***,N  . 

n  nO  n  —  nO 

Thus : 

^n(_r,t)  =  0  on  9S  ;  n  =  1,***,N 

£ 

It  has  been  shown  that  the  eigenfunctions  y^  of  the  follow¬ 
ing  problem  form  a  complete  orthonormal  set  over  the  region  S  ,  provided 
that  9S  is  a  smooth  surface: 

V  ^(r)  =  -k2  y  (r)  ;  =  0  on  9S  (2.18) 

We  should  note  that  k  =  0  is  not  an  allowed  eigenvalue  for  this  problem. 
To  show  this  let  us  assume  it  is  an  eigenvalue.  Then  there  is  an  eigen¬ 
function  y^  satisfying  the  following: 

V2yQ  =  0  ;  yQ  =  0  on  9S  .  (2.19) 

By  definition  an  eigenfunction  is  not  identically  equal  to  zero.  However 
the  function  y  =  0  does  satisfy  (2.19)  and  by  uniqueness  it  is  the  only 
solution.  Thus  there  can  be  no  eigenfunction  y^  and  zero  is  not  an 
eigenvalue.  We  also  see  from  (2.18)  that  no  y^ Or)  can  be  identically  a 
constant,  because  it  would  not  satisfy  the  boundary  condition  on  9S  . 

We  will  now  expand  each  perturbation  ij^(jr,  t)  in  a  series  of 
the  y^Or)  .  We  write 

^  (r,t)  =  l  b  (t)  y  (r)  ;  n  =  1,*",N  .  (2.20) 

n 


f 


-  12  - 


The  summation  is  over  the  allowed  eigenvalues  corresponding  to  equation 
(2.18)  in  the  region  S  . 


Now  we  will  substitute  equation  (2.20)  into  equation  (2.17) 
and  obtain: 


3b  (t)  N  df 

I  Vx)-^=  l  l  Ur 

l,  J#0  *,,(#0  j=l  \dXj 


|u^(r)bJ£(t) 


X1  X10’  " 


-  D  I  0-  b  (t)  u  (r) 
n  nX  * 


(2.21) 


n  =  !,*•*, N  . 


Now  we  multiply  both  sides  of  equation  (2.21)  by  y^(_r)  and  then 
integrate  each  side  over  the  region  S  .  Using  the  orthonormality  of 
the  y^(r_)  we  obtain: 


3b  (t)  N 
nk  _  y 

3t  ”  ^ 

3=1 


3  f 


n 


3X. 

3 


VX10>” 


b.,  (t)  -  D  kb  (t) 
j  k  n  nk 


N 

l 

3=1 


-  D  k2  6 

b  , 

Lrh 

n  nj 

3k 

X  =x  • •  •  , 

1  A10» 

J 

(2.23) 


n  =  1, 


,N 


Equation  (2.23)  is  now  a  linear  system  of  differential  equations  for  the 
b  (t)  .  It  is  similar  in  form  to  equation  (2.4).  In  analogy  with  (2.7) 
we  can  define  a  matrix  A1 (k)  and  write  this  expression  for  the  nj th 


■ 


r 


1 


13  - 


element 


9f 


n 


Anj(k)  =  3x. 


X1=X10’"‘,XN=XN0 


D  k  6 
n  nj 


(2.24) 


Then,  similar  to  equation  (2.9),  we  can  write: 


bnk^> 


N 

l 

j=i 


c  .  e 
nj 


w'(k)t 


n  =  1, * • • ,N 


(2.25) 


where  the  c^V  are  constant  amplitudes  and  the  w^  are  the  eigenvalues 
of  the  matrix  A' (k)  .  So  now  we  may  write  the  fluctuation  as: 


V-’C)  =  I  bnk(t)  Uk(-) 

n  MO.k  n  k 

N  w!  (k) t 

=  l  I  cni  e  J  Vr> 

k^0,k  j=l  J 


(2.26) 


n  =  1 ,  • • • ,N 

Again,  it  is  evident  that  whether  the  fluctuations  will  grow  in  time  or 
die  out  depends  on  whether  any  of  the  eigenvalues  of  A* (k)  possess 
positive  real  parts  or  not. 

For  the  case  of  the  boundary  conditions  given  by  (2.13)  we  may 
show  that  a  similar  analysis  will  go  through.  We  require  that  the  solu¬ 
tion  given  by  (2.15)  satisfy  the  boundary  condition.  That  is: 


dX 


§7  =  0  -  i;  (XnO  +  °"  3s  ;  n=1’"-'N 


(2.27) 


. 
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Thus  : 


difj 

(£,t)  =  0  on  3S  ,  n  =!,•••, N  .  (2.28) 


In  analogy  with  problem  (2.18)  it  has  been  shown  that  the 
eigenfunctions  of  the  following  problem  form  a  complete  orthonormal 

set  in  the  region  S  . 

y  y  3  V 

VZVk(r)  =  -kZ  Vk(r)  ;  =  0  on  9S  .  (2.29) 

It  should  be  noted  that  in  this  case  zero  is  an  eigenvalue  with  eigen¬ 
function  Vq  identically  equal  to  a  constant.  As  before,  we  expand  each 
perturbation  in  a  series  of  the  Vk  : 

^n(r>t)  =  l  ^(t)  \(£)  •  (2.30) 

Now  we  can  continue  exactly  as  before  and  we  will  obtain  the  following 
result : 

_  N  w'(t) 

=  I  l  cV  e  V  (r)  .  (2.31) 

k  j=l 

From  equations  (2.26)  and  (2.31)  we  see  that  we  obtain  very 
similar  results  for  both  types  of  boundary  conditions.  However,  in  (2.26) 
(corresponding  to  fixed  boundary  conditions)  the  summation  does  not  include 
k  =  0  .  Also  no  y  (r)  can  be  identically  equal  to  a  constant  which  means 
that  homogeneous  perturbations  are  forbidden.  On  the  other  hand,  in  (2.31) 
(corresponding  to  derivative  boundary  conditions)  the  zero  eigenvalue  is 
allowed  and  the  corresponding  eigenfunction  is  a  constant.  Thus  homogeneous 


\ 


jh  'v-' 


■ 
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perturbations  are  allowed.  The  importance  of  this  difference  will  be 
made  clearer  later  on  in  this  paper. 

The  above  arguments  are  valid  for  any  dimension  and  any  shape 
of  system,  provided  that  it  has  smooth  boundaries.  The  matrix  A'  will 
be  the  same  in  all  cases.  For  simplicity,  we  shall  work  with  one  dimen¬ 
sional  systems  in  the  rest  of  this  paper.  In  this  case  the  functions  y^ 
and  would  correspond  to  sin  kr  and  cos  kr  respectively,  where 

k  =  ~  ,  L  is  the  length  of  the  interval,  the  Z  are  integers,  including 
zero  for  cos  kr  . 

Thus  we  see  from  equation  (2.15)  that  if  the  homogeneous  steady- 
state  solution  is  unstable  and  fluctuations  do  not  regress  in  time,  then 
the  system  may  evolve  into  a  solution  that  depends  on  the  space  parameter 
_r  .  This  would  mean  that  the  concentrations  of  the  intermediate  species 
would  not  be  homogeneous  in  space.  So  this  is  one  way  that  a  chemical 
system  could  exhibit  an  inhomogeneous  state. 

It  is  clear  that  we  would  be  greatly  aided  if  we  had  a  relatively 

simple  method  for  determining  if  any  of  the  eigenvalues  of  a  matrix  have  a 

positive  real  part  or  not.  Fortunately  there  exists  a  criterion  known  as 

4 

the  "Routh-Hurwitz"  conditions  which  do  this.  We  will  now  briefly  explain 
them. 


To  find  the  eigenvalues  of  the  matrix  A  we  need  to  consider  the 


following  equation: 


. 

r 

I 


-  16  - 


'-hi 

-A12 

^3 

-A21 

w_A22 

-A23 

"A31 

"A32 

w-a33 

(2.31) 


For  an  N><N  matrix  A  this  can  be  rewritten  as  a  polynomial  of  degree 
N  or  a  "dispersion  equation": 


N  N-l  ,  N-2 

w  4-  a^w  +  a^w 


4- 


+  Vl”  +  aN  =  °  ' 


(2.32) 


The  Routh-Hurwitz  conditions  are  inequalities  among  the  coefficients 
a^,,#,,a  of  the  polynomial.  It  is  necessary  and  sufficient  that  they 
all  be  satisfied  in  order  that  the  real  part  of  every  eigenvalue  of  A 
be  negative.  Because  these  conditions  are  necessary,  it  follows  immedi¬ 
ately  that  if  one  or  more  of  them  is  violated,  then  the  real  part  of  at 
least  one  eigenvalue  is  positive.  For  future  use  the  Routh-Hurwitz 

t 

conditions  for  polynomials  of  degree  two,  three,  and  four  are  given  below 


Polynomial 


4 

w 


w  4-  a^w  4- 
3  2 

w  +  a^w  4-  a^w  4-  a^ 

3  2 

4-  a^w  +  a£W  4-  a^w  4-  a 


4 


Routh-Hurwitz  Conditions 


a2 

>  0 

>  >  0 

t 

o 

A 

a. 

>  0  ,  a.,  a~  >  a„ 

1 

’12  3 

a,  >  0  ,  an  >  0  ,  a1  >  0 
4  2  1 

2 

a3^ala2_a3^  >  ala4 


(2.33) 


\ 
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CHAPTER  3 

SPECIFIC  EXAMPLES 


Part  I 


To  begin  this  section  we  will  do  a  simple  example  to  illustrate 
the  theory  developed  in  the  previous  chapter.  We  will  consider  the  follow¬ 
ing  mechanism^: 


A  — >  X 
k2 

2X  +  Y  — >  3X 
k3 

B  +  X  — >  Y  +  D 

k4 

X  — >  E 


(3.1) 


The  overall  reaction  is  A  +  B  — >  D  +  E  .  We  have  set  the  rate  constants 
for  all  back  reactions  equal  to  zero,  so  that  the  system  is  infinitely  far 
from  equilibrium.  It  would  be  possible  to  allow  the  back  reactions^, 
although  it  would  complicate  the  algebra  substantially.  It  should  also  be 
noted  that  this  mechanism  contains  a  trimolecular  step.  While  not  realistic, 
it  is  again  simpler. 

In  writing  out  the  chemical  rate  equations  we  shall  at  first 
neglect  diffusion  and  work  only  with  the  possibility  of  a  temporal  instabil¬ 
ity.  The  chemical  rate  equations  are: 


■ 


1 


\ 
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dX 

dt 


=  k^A  4- 


k2X2Y 


-  k„BX  -  k.X 
3  4 


dY 

dt 


k3BX  -  k2X  Y 


(3.2) 


There  is  one  steady-state  solution  which  is  obtained  by 


solving : 


k,A  +  k_X  Y  -  k„BX  -  k.X  =  0 
12  3  4 


k3BX  -  k2X  Y  =  0 


(3.3) 


This  solution  is  given  by: 


xo  - 


k^A 


_  k3k4B 
Y0  =  klk2A 


(3.4) 


Now  we  shall  investigate  the  stability  of  the  steady-state 
solution.  We  shall  use  equation  (2.7)  to  find  the  matrix  A  .  From  the 
chemical  rate  equations  we  have  that: 

f,  =  k,A  +  k„X2Y  -  k„BX  -  k.X 
112  3  4 

f2  =  k3BX  -  k2X2Y 


(3.5) 


\ 
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So  from  equation  (2.7)  we  obtain: 


X0*Y0 


2k2X0Y0 


(3.6) 


The  other  matrix  elements  are  found  in  a  similar  manner. 
(3.4)  we  have  that: 


k^A 


YV 

kJ 


Using  equation 


\ 

(3.7) 

/ 


We  now  need  to  find  the  eigenvalues  of  A  .  From  equation  (2.31)  we  can 
immediately  write  this  equation  for  the  eigenvalues  of  A  : 


w  -  k_B  +  k. 
3  4 


-k 


'k  Al  2 


2  k , 


k3B 


w+k 


[k  Al2 


2  k, 


0 


(3.8) 


The  resulting  dispersion  equation  is: 


k„ 

k^A 

2 

-+*  k f  —  k^B 

+  k.k_  1 

fkiA\ 

2 

k4j 

4  3 

4  2  | 

[W 

(3.9) 


We  may  now  use  the  Routh-Hurwitz  conditions  (2.33)  for  a  second  degree 
polynomial  to  check  if  it  is  possible  to  have  an  instability.  We  have: 


20  - 


a2  k4k2 


k  A 


(3.10a) 


al  =  k2 


lc  A.2 

—  +  k4  -  k3® 


(3.10b) 


It  is  obvious  that  a^  >  0  for  all  physically  meaningful  values  of  the 

rate  constants  and  A  .  However  a^  is  not  greater  than  zero  for  all 

reasonable  values.  We  see  that  there  exists  a  critical  value  of  B  ,  call 

it  B  ,  such  that  for  B  >  B  ,  a.  <  0  .  This  critical  value  is  given 
c  c  l 

by : 


B 

c 


(3.11) 


Thus  for  B  >  B^  ,  the  steady-state  solution  to  the  chemical  rate  equations 
(3.2)  becomes  unstable.  It  can  be  expected  that  if  it  is  possible  the 
system  will  evolve  to  a  time  dependent  solution  for  B  >  Bc  . 


It  may  be  noted  that  the  chemical  rate  equations  (3.2)  are  non¬ 
linear.  This  will,  in  general,  be  the  case  if  we  wish  to  find  a  meaningful 
instability  in  the  chemical  rate  equation.  A  simple  example  will  show  why 
linear  equations  would  be  unsatisfactory.  Consider  the  following  mechanism 


A  — >  X 
k2 

B  +  X  —  >  2X 
k3 

X  — >  C 
k4 

D  +  X  — >  E 


(3.12) 


4  I 
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It  is  convenient  to  set  all  rate  constants  equal  to  one.  The  chemical 
rate  equation  is: 


dX 

dt 


=  A  +  (B-l-D)X 


(3.13) 


We  have 


X0  B-l-D 


(3.14) 


The  matrix  A  in  this  case  consists  of  the  single  element  (B-l-D)  . 

The  condition  for  instability  is  B  <  1  +  D  .  However,  we  see  from 
equation  (3.14)  that  when  the  instability  condition  holds,  the  steady- 
state  solution  is  negative  which  does  not  make  sense  physically.  Similar 
phenomena  appear  to  happen  in  systems  which  contain  more  than  one  inter¬ 
mediate  species  and  have  linear  chemical  rate  equations. 


Now  we  shall  demonstrate  the  existence  of  an  instability  that 
results  from  fluctuations  which  depend  on  both  time  and  space  parameters. 
We  again  consider  mechanism  (3.1)  but  this  time  in  the  chemical  rate 
equations  we  allow  for  diffusion^  *  ye  have: 


dX 

dt 


=  k^A  + 


k2X2Y 


-  k~BX  -  k.X  +  D 
3  4  x 


3Y 

9t 


k„BX  -  k0X2Y  +  D 
3  2  y  .  2 

dr 


(3.15) 


\ 

' 

( 


■ 
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The  steady-state  homogeneous  solution  is  still  given  by 
equations  (3.4).  This  time,  however,  in  considering  the  stability  of 
this  solution  we  must  use  the  matrix  A' (k)  whose  elements  are  given 
by  equation  (2.24).  f  and  f ^  are  still  given  by  equations  (3.5)  and 

(3.6)  respectively.  We  easily  obtain  for  A* (k) : 


The  corresponding  equation  for  w  is 


w  -  k~B  +  k.  +  k  D 
3  4  x 


-k. 


k^A 


k3B 


w  +  k 


fkiA 

2  l  k. 


2 

+  k  D 


=  0  (3.17) 


This  may  be  rewritten  as: 


2 

w 


4-  w(k^ 

+  k4k2 

+  k2k, 
4 


+  k. 


[k  A] 


-  k„B  +  k2  D  +  k2  D  ) 
3  x  y 


k1A,2 


+  k  D  k_ 
x  2 


k1A»2 


k  k„B  D 

3  y 


D  +  k4  D  D  = 
Y  x  y 


0 


(3.18) 
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For  algebraic  simplicity  we  shall  at  this  point  set  all  rate  constants 
equal  to  one.  Then  equation  (3.18)  becomes: 


2  2  2  2  2 

w  +  w  (1  +  A  -  B  +  k  D  +  k  D  ) 

x  y 


+  A2  +  k2A2D  -  k2  BD  +  k2D  +  k4D  D  =  0 
x  y  y  •  x  y 


(3.19) 


The  Routh-Hurwitz  stability  condition  of  a^  >  0  gives  us  the  following 
critical  value  of  B  : 


B  =  1  +  A2  +  k2(D  +D  ) 
c  x  y 


(3.20) 


When  B  >  Bc  ,  the  homogeneous  steady-state  solution  will  become  unstable. 
However,  it  is  now  possible  to  violate  the  condition  that  a?  >  0  .  This 
gives  us  another  critical  B  which  we  will  denote  by  B^  : 


b;oo  = 


A2  D  A  , 

+  1  +  k2D 

2_  D  x 


(3.21) 


k  D 


It  is  of  interest  to  find  which  value  of  k  gives  the  smallest 

2 

B'  .  To  find  this  we  differentiate  B*  with  respect  to  k  and  set  the 
C  c 

resulting  expression  equal  to  zero.  We  obtain  the  following  equation: 


j-  ~  k4  D  =0 
D  x 

y 


(3.22) 


This  is  easily  solved  to  yield: 


k2  =  _ A _ 

C  (D  D  )1/2 
x  y 


(3.23) 


\ 


. 


■ 
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Putting  this  value  into  equation  (3.21)  yields: 


B' 

c 


(3.24) 


At  this  point,  Prigogine  who  developed  this  theory  concludes 
by  stating  that  for  B  greater  than  the  minimum  of  B^  and  B^  ,  the 
homogeneous  steady-state  solution  becomes  unstable.  The  critical  value 
of  B  corresponds  to  a  "symmetry  breaking  transition".  If  beyond  the 
instability,  a  space  dependent  state  exists,  Prigogine  calls  it  a  "dissi¬ 
pative  structure". 

However,  it  is  interesting  to  examine  more  closely  the  role 
the  boundary  conditions  may  play.  First  let  us  consider  the  derivative 
boundary  conditions.  For  this  case  we  have: 


3X(r,t) 


3r 


 9X(r, t) 


r=0 


3r 


=  0 


r=L 


(3.25) 


3Y(r>t) 

3r 


r=0 


9Y(r,t) 


3r 


r=L 


0  . 


With  this  type  of  boundary  condition  it  was  shown  previously 
that  it  is  possible  to  have  a  perturbation  independent  of  r  .  That  is, 
we  may  have : 


X  =  XQ  +  \p(  t) 


Y  =  YQ  +  4>(t) 


(3.26) 


- 


f 
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This  homogeneous  perturbation  corresponds,  of  course,  to  k  =  0  .  At 
k  =  0  ,  Bc  takes  on  its  minimum  value  which  is: 


B  (min)  =  1  +  A 
c 


(3.27) 


If  the  values  of  D  and  D  are  such  that: 

x  y 


B  (min)  =  1  +  A 
c 


B'  = 
c 


(3.28) 


then  the  instability  will  always  be  caused  first  by  a  homogeneous  pertur¬ 
bation  such  as  in  (3.26).  Initially  the  system  will  go  into  a  time 
dependent,  space  independent  solution. 


On  the  other  hand,  if  the  boundary  conditions  are  of  the  other 
type,  that  is: 


X(r,t) 


r=0 


=  X(r,t) 


r=L 


=  X 


0 


Y(r,t) 


r=0 


=  Y(r,t) 


r=L 


=  Y 


0 


(3.29) 


a  different  situation  exists.  For  these  boundary  conditions,  perturbations 
of  the  form  of  (3.26)  are  impossible  to  have.  Since  we  cannot  have  pertur¬ 
bations  which  correspond  to  k  =  0  ,  B^  will  take  on  its  minimum  value 
2 

when  k  is  the  smallest  possible.  Since  the  eigenfunctions  ]i^  are 
sin  kr  ,  this  means  when  k  =  —  . 


In  this  case  we  have  that : 


\ 
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B=l  +  A2+(D+D)^r 
C  X  y  L2 


(3.30) 


The  instability  will  be  caused  by  whichever  is  the  smaller: 


1  +  A2  +  (D  +  D  )  ~  or  B' 

x  y  l2  c 


[1  +  (irfAi2 


(3.31) 


However,  no  matter  whether  B  or  B1  is  the  smaller,  the  fluctuation 

c  c 

must  be  inhomogeneous  and  beyond  the  instability  the  system  evolves  into 
a  space  dependent  state. 

Thus  we  see  that  there  is  a  fundamental  difference  between  the 
two  types  of  boundary  conditions.  With  one  of  them,  homogeneous  pertur¬ 
bations  are  allowed  and  with  the  other  homogeneous  perturbations  are  not 
allowed.  The  importance  of  this  difference  is  made  clearer  by  a  consider¬ 
ation  of  the  following  mechanism.  Prigogine^  has  stated,  without  proof, 
that  this  mechanism  may  exhibit  a  symmetry  breaking  transition.  The 
mechanism  is : 


A  +  X  — >  2X 
k2 

X  +  Y  —  >  2Y 
k3 

Y  -I-  V  —  >  V' 


(3. 32) 


k, 

V'  — >  E  +  V 


The  overall  reaction  is  A  E  .  The  chemical  rate  equations  for  this 
system  are: 


\ 
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3X 

3t 


3y 

at 


3V 

3t 


3V' 

at 


a2x 

=  AX  -  k„XY  +  D  — £ 

2  x  ~  2 

3r 

a2Y 

=  k0XY  -  k_YV  +  D  — i- 
2  3  y  ar2 

2 

=  k.V'  -  k  YV  +  D 

4  3  v  „  2 


k_YV  -  k.V'  +  D  , 
3  4  v 


2 

a  V' 

a  2 
3r 


(3.33a) 

(3.33b) 

(3.33c) 

(3.33d) 


Since  the  product  k^A  always  occurs  together,  it  is  convenient 
to  denote  it  simply  by  A  .  This  convention  will  always  be  followed  where 
applicable  in  the  remainder  of  this  paper. 


The  homogeneous  steady-state  solution  to  the  system  of  equations 


(3.33)  is 


k  V' 
K4V0 


L0 


Y„  = 


0  k. 


(3.34) 


V0  = 


k4k2V0 

k3A 


Vq  is  an  arbitrary  parameter  which  arises  because  the  set  of  algebraic 
equations  of  which  (3.34)  is  the  solution  has  two  linearly  dependent 
equations.  (See  equations  (3.33c)  and  (3.33d)).  We  wish  to  test  the 
stability  of  this  homogeneous  steady-state  solution.  Using  equations 
(2.24)  and  (2.31)  we  easily  obtain  the  following  equation: 


. 
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w  +  kzD 


x 


■k.  A 


0 


0 


WS 

A 


w  +  k  D 


k2k4V0 


w  + 


0 


V 


k3A  2 
-=-  +  k  D 
k^  v 


-k3A 


0 


0 


-k, 


w  +  k.  +  k2D' 
4  v 


=  0  (3.35) 


It  is  easier  to  see  what  is  occurring  if  we  denote  groups  of  constants  in 
the  following  way: 


k3A 


=  K 


k2k4V0 


B 


k2D 


x 


=  a 


;  k2D  =  b 

y 

,  2_ 

;  k  D  =  c 

v 

;  k2D'  =  d 
*  v 


(3.36) 


Now  equation  (3.35)  can  be  written  as: 


w  +  a 


-A 


B 


w  +  b 


0 


K 


0 


0 


B 


-B 


w  +  K  +  c 


-K 


0 


0 


-k. 


w  +  k.  +  d 
4 


=  0 


(3.37) 


\ 
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Expanding  and  simplifying  gives  us: 

4  3  2 

w  +  w(K+c  +  k^-fd  +  b+a)  +  w  (aK  +  ac  +  ak^  +  ad  +  Kd) 

2 

+  w  (ck^  +  cd  4-  AB  +  bK  +  be  +  bk^  4-  bd  —  KB) 

4-  w(aKd  +  ack.  4-  acd  +  abK  +  abc  4-  abk.  4-  abd 
4  4 

+  ABK  +  ABc  +  ABk.  +  ABd  +  BKd  +  bck.  +  bed  -  aKB) 

4  4 

+  abKd  +  abek.  +  abed  -  aKBd  +  ABKd  +  ABek.  +  ABed  =  0  .  (3.38) 

4  4 

Let  us  consider  first  the  derivative  boundary  conditions.  In 
this  case,  as  shown  above,  we  must  allow  for  homogeneous  perturbations. 
That  is; 


a=b=c=d=0 


(3. 39) 


In  this  case  equation  (3.42)  simplifies  considerably.  It  becomes: 


/  o  o 

w  +  w  (K+k4)  +  w  (AB-KB)  +  w(ABk4  +  ABK)  =  0 


(3.40) 


We  immediately  see  that  one  solution  is  w  =  0  .  Factoring  this  out 


leaves  us  with: 


3 

w 


+  w 


(K+k4)  +  w (AB-KB)  +  ABk4 


+  ABK  =  0 


(3.41) 


Now  we  apply  the  Routh-Hurwitz  condition  for  stability  to  equation  (3.41). 
The  condition  that  a^a2  >  a3  becomes: 


y 
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(K+k4) (AB-KB)  >  AB(K+k4) 


(3.42) 


or : 


A  -  K  >  A  . 


(3.43) 


In  terms  of  the  natural  parameters  of  the  system  this  is: 


A  -  >  A  or 


(3.44) 


Quite  obviously,  this  can  never  be  satisfied  since  k^  and  k4  must 


both  be  positive.  This  means  that  with  respect  to  homogeneous  perturba¬ 
tions  the  homogeneous  steady-state  solution  is  always  unstable  for  all 


values  of  the  rate  constants  and  concentrations  of  A  and  .  Clearly 


a  symmetry  breaking  transition  is  not  possible  under  the  derivative  boundary 
conditions,  since  the  steady-state  homogeneous  solution  is  never  stable. 


Now  we  will  consider  for  fixed  boundary  conditions  the  stability 


of  the  steady-state  homogeneous  solution  given  by  equations  (3.34).  In 
this  case  it  is  easy  to  show  that  there  do  exist  values  of  the  rate  constants, 
concentrations,  diffusion  constants,  and  box  size  such  that  the  homogeneous 
steady-state  solution  is  stable  to  all  allowed  fluctuations.  As  an  example. 


choose  the  rate  constants,  A  and  such  that  K  =  1  and  B  =  1  .  In 


addition  choose  all  of  the  diffusion  constants  equal  to  one  and  choose  the 
size  of  the  box  to  be  TT  .  This  makes  the  smallest  allowed  k  value  equal 
to  one.  Then  it  is  quite  easy  to  show  that  all  of  the  Routh-Hurwitz  condi¬ 
tions  for  stability  will  hold  in  this  case  for  any  allowed  fluctuations.  On 
the  other  hand,  it  is  also  possible  to  show  that  for  other  values  of  the 


. 
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parameters  a  symmetry  breaking  transition  is  possible.  Let  us  consider 
the  Routh-Hurwitz  condition  a^  >  0  and  show  that  this  may  be  violated 
under  certain  conditions.  For  convenience,  we  will  set  all  rate  constants 
equal  to  one.  Then  using  (3.38)  we  have: 

V' 

K  =  A  ;  B  =  -£■  (3.45) 

r 

a,  =  k  (D  D  D'A  +  D  D  D  +  D  D  D  D'k  ) 

4  xyv  xyv  xyvv 

+  V'  k2(AD’  +  D  +  k2D' (D  -  D  ))  .  (3.46) 

0  v  v  v  v  x 


Equations  (3.36)  have  also  been  used  to  obtain  the  above  equations.  If 
a^  <  0  we  know  that  the  homogeneous  steady-state  solution  is  unstable. 
That  is,  if: 


< 


/  Q 

-k  (D  DD'Dk  +DDD'A+DDD) 
_ xyvv _ xyv _ xyv 

[AD'  +  D  +  k2D'  (D  -  D  )] 

V  V  V  V  X 


(3.47) 


for  the  denominator  positive.  If  the  denominator  is  negative  (which 
implies  >  D  )  we  have  the  following  condition  for  instability: 


V'  > 
0 


/  o 

-k  (D  DD'Dk  +DDD'A+DDD) 
xyvv _ xyv _ xyv 

[AD'  +  D  +  k2D' (D  -  D  )] 

V  V  V  V  X 


=  F(k) 


(3.48) 


The  right  hand  side  of  (3.47)  is  negative,  so  it  would  be  impossible  for 
Vq  to  satisfy  this.  However  for  certain  kjV^  can  satisfy  (3.48).  It 
is  instructive  to  consider  a  rough  sketch  of  F(k)  .  (See  Figure  1.)  The 
region  where  k  <  k^  corresponds  to  equation  (3.47).  The  value  of  k^ 


" 
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can  be  found  by  differentiating  F(k)  with  respect  to  k  and  setting 
the  resulting  expression  equal  to  zero.  Solving  this  equation  will  give 
k  .  If  we  substitute  for  k  by  k  into  F(k)  given  by  (3.48)  we 
have  the  minimum  value  that  must  be  greater  than  in  order  that  the 

steady-state  homogeneous  equation  be  unstable.  That  is,  instability  occurs 
if: 


Vq  >  F(kc)  .  (3.49) 

Thus  for  boundary  conditions  that  fix  the  concentrations  of  X,Y,V,V’  on 
the  boundary  it  is  possible  to  have  a  symmetry  breaking  transition. 

In  summary  we  see  that  whether  or  not  a  system  can  have  a  symmetry 
breaking  transition  may  be  determined  by  the  type  of  boundary  conditions 
imposed  on  it.  The  mechanism  we  have  just  investigated  exhibited  this 
phenomenon  for  fixed  concentrations  at  the  boundary,  but  did  not  when  the 
boundary  conditions  required  only  that  the  normal  derivative  of  the  concen¬ 
trations  vanish  at  the  boundary. 


f 
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Part  II 


The  chemical  rate  equations  are,  in  general,  a  set  of  coupled, 
non-linear,  partial  differential  equations.  It  would  usually  be  diffi¬ 
cult  to  study  the  actual  time  dependent,  space  dependent  solution  to  this 
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set  of  equations.  For  this  reason,  Lefever  *  *  has  devised  a  model 

which  simplifies  the  analysis.  We  will  examine  this  model  for  the  mechan¬ 
ism  given  by  (3.1).  We  assume  that  the  reaction  takes  place  in  two  identi¬ 
cal  boxes  labelled  one  and  two.  The  concentrations  of  A  ,  B  ,  D  ,  and  E 
are  fixed  and  are  the  same  in  both  boxes.  The  concentrations  of  X  in  box 
1,  denoted  by  ,  is  constant  but  need  not  equal  X2  .  A  similar  situation 

is  true  for  and  Y2  •  X  and  Y  are  allowed  to  diffuse  between  the  two 

boxes . 


The  chemical  rate  equations  for  this  model  are: 


dX 


—  =  A  +  k0X^  Y.  -  BX.  -  k, X  +  D  (X9  -  X..  ) 
dt  211  1  41  x  2  1 


dX 


~  -  A  +  k9X9  Y_  -  BX0  -  k.X0  +  D  (X  -  X9) 

dt  222  2  42  xl  2 


(3.50) 


dY. 


"dt"  “  BX1  k2Xl  Y1  +  VY2 


V 


dY, 


"dt  =  BX2 


k2X2  Y,  +  D  (Yl  -  Y2) 


Note  that  equations  (3.50)  are  ordinary  differential  equations,  not 
partial.  This  is  why  this  model  is  simpler  to  consider. 


\ 
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The  homogeneous  steady-state  solution  is: 

X10  ■  X20 
Y10  =  Y20 

The  inhomogeneous  steady-state  solution  is  any  X^  ,  X2q  ,  Y^q  ,  Y2q 
which  satisfy  the  following  algebraic  equations: 


k.B 

4 

k2A 


(3.51) 


A  +  k2  X10  Y10  BX10  k4X10  +  Dx(X20  X10^  ° 


A  +  k2  X20  Y20  '  BX20  -  k4X20  +  Dx(X10  '  X20}  =  0 


BX10  k2  X10  Y10  +  °y(Y20  Y10}  ° 


BX20  k2  X20  Y20  +  Dy(Y10  Y20')  ° 


(3.52) 


and 


L10 


¥  X 


20 


10 


¥  Y 


20 


In  order  to  be  able  to  study  the  stability  of  both  homogeneous 
and  inhomogeneous  solutions  to  equations  (3.50),  Lefever  chose  all  rate 
constants  equal  to  one  and: 

D  =  1  ;  A  =  2  .  (3.53) 

X 

Doing  this  will  leave  us  with  only  D  and  B  as  arbitrary  parameters. 


\ 


■ 


36  - 


Let  us  consider  first  the  stability  of  the  homogeneous  steady- 
state  solution.  We  could  use  equations  (2.7)  and  (2.31)  to  give  us  a 
4x4  matrix  A  and  thus  a  fourth  order  dispersion  equation  to  analyze. 
However,  it  is  possible  to  consider  instead  two  matrices,  each  of  which 
is  2x2  .  First  we  will  consider  only  homogeneous  perturbations  to  the 
homogeneous  solution.  That  is: 


Xi  =  Xi0  +  ^(t)  ;  1  =  1»2 
Yi  =  Yi0  +  *(t>  ;  1  *  1>2 


(3.54) 


Solution  (3.54)  is  slightly  different  from  the  perturbed  solution  considered 
in  chapter  2.  Here  the  fluctuations  are  not  completely  general.  Instead 
we  have  set : 


^(t)  =  ip2(t)  =  ijj(t)  and  4>1(t)  =  <f>2(  t)  =  cj)(t) 


Thus  we  have  X?  -  =  0  and  Y?  -  =  0  .  Substituting  equations 

(3.54)  into  (3.50)  and  linearizing  gives  the  following  set  of  equations: 


dip(t) 

dt 


(B-l)  ip(t)  +  4<|>(t) 


(3.55) 


=  BtKt)  -  4(f)  (t) 


(3.56) 


We  have  used  here  the  fact  that  all  rate  constants  are  one  and  equations 
(3.51)  and  (3.53).  The  matrix  A  corresponding  to  this  system  is: 


. 
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A  = 


B  -  1 


\ 


-B 


(3.57) 


So  the  equation  for  the  eigenvalues  of  A  is: 


w  -  B  +  1 


+B 


-4 


w  +  4 


=  0 


(3.58) 


The  resulting  dispersion  equation  is: 


w  +  w(5-B)  +4=0 


(3.59) 


Using  the  Routh-Hurwitz  conditions  we  find  that  the  homogeneous  steady- 
state  solution  is  unstable  to  homogeneous  perturbations  if: 


B  >  5 


(3.60) 


Now  let  us  consider  the  stability  of  the  homogeneous  solution 
with  respect  to  inhomogeneous  perturbations.  We  write  for  the  perturbed 
solutions : 

xi =  xio +  h(t) 

X2  -  X20  +  *2(t) 

Y1  =  Y10  +  VC) 

Y2  =  Y20  +  *2(t)  ' 


(3.61) 


/ 
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For  inhomogeneous  perturbations  we  require  ^ (t )  ^  ^(t)  and 
^(t)  ?  <|>2(t)  ’  nOW  su^st;*-tute  (3-61)  into  (3.50),  linearize,  and 

use  equations  (3.51)  and  (3.53)  as  well  as  the  fact  that  the  rate 
constants  are  one.  This  yields: 


d^1 

“dt 


=  (B-l)^  +  4<f>^  +  (\p2  -  ip^ 


\ 


dip2 

~dt 


=  (B-l)ip2  +  4<J>2  -  (\p2  -  \p±) 


dh  _  _  _  x 

-at  -  B*i  -  +  y*2  -  V 

d^2  -  -  -  - 

—  -  B^2  -  4^  -  Dy(^  - 


(3.62) 


(3.63) 


Equations  (3.62)  can  be  combined  to  give  one  equation  for  (ip2  -  ip^)  . 
A  similar  arrangement  can  be  performed  with  equation  (3.63).  Thus  we 
obtain : 


d(ip?  ~  ip,)  _ 

- (B  -  3)  (^2  ~  V  +  *(*2  -  V 

(3.64) 

d(<P2  -  40  _  _ 

— -  -a-  - =  BOP  -  ip  )  -  (4  +  2D  )  (4>  -  4  )  . 


We  now  use  the  system  of  equations  given  by  (3.64)  to  determine  the  matrix 
A  .  It  is  easy  to  see  that  the  equation  for  the  eigenvalues  of  A  is 
given  by : 


■ 


■ 
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w  -  B  +  3 


B 


-4 


w  +  4  4-  2D 

y 


=  o 


(3.65) 


From  this  we  obtain  the  following  dispersion  equation: 

2 

w  +  w(7  -  B  +  2D  )  +  6D  4-  12  -  2BD  =  0  .  (3.66) 

y  y  y 

Using  the  Routh-Hurwitz  conditions  we  see  that  the  homogeneous  steady- 

state  solution  is  unstable  with  respect  to  inhomogeneous  perturbations  if: 

B  >  7  +  2D  (3.67a) 

y 

or 

3D  +  6 

B  >  B  =  — -  .  (3.67b) 

c  D 

y 

Which  of  these  is  reached  first  will,  of  course,  depend  on  the  value  of 
D^  .  However,  before  (3.67a)  is  satisfied  the  homogeneous  steady-state 
solution  would  become  unstable  with  respect  to  homogeneous  fluctuations. 
(See  equation  (3.60)). 

Next  we  shall  look  at  the  stability  of  inhomogeneous  steady- 
state  solutions  to  the  kinetic  equations.  In  this  case  it  is  impossible 
to  simplify  the  analysis  and  we  must  consider  a  fourth  order  dispersion 
equation.  Using  equations  (2.7)  and  (2.31)  we  obtain: 


' 


\ 


■ 


-  40  - 


-1 


0 


-1 


0 


0  (3.68) 


0 


-D 


y 


0 


-D 


y 


Here  ,  X^q  ,  Y^q  ,  Y^q  is  a  solution  to  the  set  of  equations  given 


by  (3.52).  These  equations  must  be  solved  numerically  for  given  values 
of  B  and  (using  the  values  already  specified  for  the  rate  constants, 

A  and  ).  Then  equation  (3.68)  can  be  solved  to  determine  the  stability 

of  the  solution  or  the  dispersion  equation  can  be  obtained  from  (3.68)  and 
the  Routh-Hurwitz  conditions  used  to  check  for  stability. 

Lefever  has  carried  out  the  numerical  analysis  and  his  results 
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are  shown  in  Figure  2.  We  quote  from  his  paper  to  describe  Figure  2 
"The  homogeneous  steady-state  is  stable  in  domain  I.  It  becomes  unstable 
with  respect  to  inhomogeneous  perturbations  beyond  the  curve  b  .  [It  is 
also  unstable  for  B  >  7  +  .  Lefever  apparently  did  not  consider  this 

latter  possibility  and  so  did  not  draw  in  the  appropriate  line.]  In 
regions  II,  III,  and  V  no  time-independent  stable  state  exists.  On  the 
contrary,  in  region  IV  the  inhomogeneous  steady-state  is  stable....  This 
region  corresponds  therefore  to  what  was  called  a  dissipative  structure." 

In  this  chapter  we  have  reviewed  the  works  of  Prigogine,  Lefever, 
and  Nicolis  which  show  that  the  homogeneous  steady-state  solution  to  the 
chemical  rate  equations  can  become  unstable  with  respect  to  spatial  fluctua¬ 
tions.  This  implies  that  the  stable  solution,  if  there  is  one,  may  be  space 


. 

•  < 


~4! 


FIGURE  2 
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dependent.  We  have  also  reviewed  the  work  of  Lefever  on  the  two  box 
model  which  provides  further  strong  evidence  for  the  existence  of  a 
space  structure.  In  addition,  the  importance  of  the  boundary  conditions 
was  examined,  which  had  not  been  done  previous  to  this  paper.  An  example 
was  presented  which  showed  that  the  type  of  boundary  conditions  imposed 
on  the  system  can  determine  whether  or  not  a  symmetry  breaking  transition 
is  possible. 
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CHAPTER  4 

MULTIPLE  STEADY-STATES 


So  far  the  chemical  rate  equations  have  had  only  one  steady- 
state  homogeneous  solution.  However,  since  the  equations  are  non-linear, 
it  is  possible  to  have  more  than  one  solution.  The  following  simple 
mechanism  exhibits  this  phenomena  and  we  will  study  it  in  detail: 


A  +  2X 


3X 


B  +  X 


C 


(4.1) 


We  wish  to  look  at  this  system  in  the  two  box  model.  Using 
this  model,  the  chemical  rate  equations  are: 


dXx 

dt 


ax;  - 


k-ixi 


BX  +  C  +  dx(x2  ~  xx) 


dX2 

“dt 


=  AX„  - 


k_iX2 


BX2  +  C  +  D  (X  -  X2) 


(4.2) 


For  the  time  being  we  will  designate  any  homogeneous  steady-state  solution 
by  Xq  .  We  will  now  study  its  stability.  Using  equations  (2.7)  and 
(2.31)  we  obtain  the  following  equation: 


■ 
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w 


+  3k_±X 


2 

0 


+  B  +  D 
x 


w  -  2AXq  + 


3k  ,Xn  +  B  +  D 
-10  x 


0 


(4.3) 


This  can  be  rewritten  as: 

w2  +  2w(B  +  D  +  3k  X2  -  2AX  ) 

X  -L  U  U 

+  (B  +  3k  ,  X2  -  2AKn  +D)2-D2=0  .  (4.4) 

-1  U  Ox  x 

Using  the  Routh-Hurwitz  conditions  we  see  that  the  conditions  for  stability 
are  : 

B  +  3k  ..X2  -  2AXn  >  -D  (4.5a) 

-1  0  0  x 

B  +  3k  X2  -  2AX_  >  -2D  (4.5b) 

-10  0  x 

B  +  Sk^X2  -  2 AXq  >  0  .  (4.5c) 


If  (4.5c)  is  satisfied  the  other  two  will  automatically  be 
satisfied.  So  equation  (4.5c)  is  the  only  equation  we  will  have  to  consider 
in  the  stability  analysis  of  a  homogeneous  steady-state  solution. 


To  find  the  homogeneous  steady-state  solution  we  need  to  solve 
the  following  equation  for  Xq  : 


X 


3 

0 


+  BX 


0 


C  =  0 


(4.6) 


i 


\ 
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Since  this  is  a  cubic  equation  there  are  three  roots  possible,  thus 
allowing  multiple  steady-states.  In  its  present  form  equation  (4.6)  is 
difficult  to  solve.  To  facilitate  solving  it  we  will  introduce  the 
parameters  a  ,  c  ,  d  ,  and  e  ,  through  the  following  implicit  equations 


k-l = 

c 

(4.7a) 

A  = 

ac  4-  d 

(4.7b) 

B  = 

e  +  ad 

(4.7c) 

C  = 

ea 

(4. 7d) 

Now  equation  (4.6,)  may  be  written  as: 

(XQ  -  a) (cXq  -  dXQ  +  e)  =  0  .  (4.8) 


The  solutions  to  this  equation  are  immediate.  We  find: 


xo  =  a 


X 


/•2 


0+ 


+d  +  /  d  -  4ce 
2c 


+d 


-  /-2 


d  -  4ce 


2c 


(4.9a) 


(4.9b) 


(4.9c) 


Depending  on  the  choice  of  parameters  c  ,  d  ,  and  e  ,  only  (4.9a)  may 
be  real  and  positive  and  thus  physically  meaningful  or  all  three  of  the 
solutions  may  be  meaningful. 


\ 
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Now  we  will  use  condition  (4.5c)  to  study  the  stability  of 
these  three  solutions.  It  is  convenient  to  use  equations  (4.7)  to 
rewrite  condition  (4.5c).  It  becomes: 


e  +  ad  +  3cXq  -  2(ac+d)  XQ  >  0 


(4.10) 


First  let  us  test  the  stability  of  the  X^  =  a  solution.  We  obtain: 


e  -  ad  +  ca  >0 


(4.11) 


This  is  a  quadratic  equation  in  a  .  The  stability  condition  (4.11)  may 
be  rewritten  as  the  following  two  equations: 


a  > 


d  + 


(4.12a) 


a  < 


-  /-2 


d  -  4ec 


2c 


=  X 


0- 


(4.12b) 


Written  this  way  we  have  the  rather  surprising  result  that  as  a  is 
varied  while  d  ,  e  and  c  are  held  constant,  the  solution  X^  =  a 
becomes  unstable  as  it  passes  through  the  values  of  the  other  two  homogen¬ 
eous  steady-state  solutions. 


Now  we  will  consider  the  stability  of  the  solution: 


d  + 


-  4ec 


2c 


This  solution  comes  from  solving: 


\ 
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+  e 


0 


(4.13) 


Rearranging  gives  us: 


e 


Substituting  this  into  condition  (4.10)  yields: 


(4.14) 


-X^+d  +  2e  -  da  +  2  X^+ca  <  0 


or 


a(2XQ+c-d)  <  X^d  -  2e 


(4.15) 


Substituting  for  X^+  gives: 


■2  •  ■  A2 


.  d  +  d  /  d  -  4ec  -  4ec 
a  <  -  =  X 


/-2 


0+ 


2c  /  d  -  4ec 


(4.16) 


In  a  similar  way  we  find  that  the  condition  for  stability  of  the  X 
solution  is  : 


0- 


a  >  XQ_  .  (4.17) 

Thus  we  see  that  when  d  ,  e  ,  and  c  are  chosen  such  that 
three  solutions  are  possible,  the  stability  of  each  solution  is  determined 
by  the  value  of  a  with  respect  to  d  ,  e  ,  c  .  It  is  quite  interesting 
that  just  as  one  solution  becomes  unstable,  another  becomes  stable.  By 
looking  at  equations  (4.12),  (4.16),  and  (4.17)  we  see  that  there  are 
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always  two  stable  solutions.  (The  question  of  stability  cannot  be 
discussed  if  a  =  or  if  a  =  (See  equations  4.12.)  because  then 

the  theorem  that  assured  us  our  linearization  process  was  valid  no  longer 
applies. ) 

The  solution  X^  =  a  is  of  particular  importance  because  it  is 
the  solution  at  thermodynamic  equilibrium.  At  equilibrium  the  rates  of 
forward  and  backward  reactions  are  equal.  That  is: 


2  3 

A  X  =  k  .X 
e  -1  e 


(4.18a) 


and 


B  X  =  C 
e 


(4.18b) 


Thus  : 


Using  equations  (4.7)  we  have: 


X 

e 


ea  _  ac  +  d 
e  +  ad  c 


This  can  be  satisfied  only  if: 


d  =  0 

X  =  a 
e 


(4.19) 


(4.20) 


(4.21a) 

(4.21b) 


Equation  (4.21a)  assures  us  that  there  is  only  one  physically  meaningful 


i 


1 


\ 
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solution  since  the  other  two  will  now  be  imaginary  quantities.  Because 
of  the  physical  significance  of  the  Xq  =  a  solution,  it  is  not  too 
surprising  that  the  parameter  a  seems  to  play  such  an  important  role  in 
the  stability  conditions. 

After  mechanism  (4.1)  had  shown  that  there  could  exist  multiple 

steady-state  homogeneous  solutions  to  kinetic  equations,  another  example 

22 

done  by  Edelstein  appeared  in  the  literature  .  This  example  also  intro¬ 
duced  thenultiple  solutions  through  a  cubic  equation,  but  the  mechanism 
governing  the  system  contained  two  intermediate  species  instead  of  only  one 
as  considered  in  this  paper. 

Mechanism  (4.1)  also  allows  us  to  illustrate  another  interesting 
phenomenon.  In  the  previous  chapter  the  existence  of  a  spatially  dependent 
solution  was  shown  by  proving  that  the  homogeneous  steady-state  solution  could 
become  unstable.  However,  it  may  also  be  possible  for  a  non-homogeneous 
solution  to  be  stable  simultaneously  with  a  homogeneous  steady-state  solu¬ 
tion^.  If  this  happened,  the  actual  state  of  the  system  would  be  determined 
by  the  initial  conditions  imposed  on  the  system.  Different  initial  conditions 
could  conceivably  lead  to  different  stable  states.  We  will  show  that  with  the 
two  box  model  it  is  possible  to  have  a  homogeneous  steady-state  stable  solu¬ 
tion  and  an  inhomogeneous  steady-state  stable  solution  simultaneously. 

Because  of  the  difficulty  of  solving  analytically  for  the  inhomogeneous  steady- 
state  solution  we  will  give  a  numerical  example. 

The  algebraic  equations  for  an  inhomogeneous  steady-state  solution 


are : 


r 


M  . 
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(X10  -  a)(cXl0  -  <“10  +  e)  =  Vx20  -  X10> 

(4.22) 

(X20  ”  a)(cX20  "  dX20  +  e)  =  "Dx(X20  “  X10}  * 

First  we  must  determine  the  criterion  for  stability  of  the 
inhomogeneous  steady-state  solution.  By  analogy  with  equation  (4.3)  we 
have : 


w  -  2AX..  A  +  3k  X  +  B  +  D 
ID  -1  10  x 


-D 


-D 


w  - 


X 


2 AX  +  3k  -X-.  +  B  +  D 
20  -1  20  x 


=  0 


Expanding  gives : 


2  2 

w  +  w(3k  1X_n+2B+2D  -2AX. n-2AX. A) 
-1  10  x  10  20 


(4.23) 


+  (3k  _xJn+B+D  -2AX  n)  (3k  ,+B+D  -2A XnJ  -  D  2  =  0 

-1  10  x  10  -1  x  20  x 


(4.24) 


Using  the  Routh-Hurwitz  conditions  and  equation  (4.7)  we  can  write  the 
conditions  for  stability  as: 


3c(X20+X2q)  +  2e  +  2ad  -  2  (ac+d)  (X10+X2())  +  2^  >  0 

(3cX2  ■ Fe+ad+D  -2 (ac+d)X1 A) (3cX2  +e+ad+D  -2(ac+d)XOA)  -  D  2  >  0 
10  x  10  20  x  20  x 


(4.25a) 

(4.25b) 


We  will  choose  the  following  values  for  our  parameters: 
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D 


1 


a  =  2 


x 


(4.26) 


d  =  +4c 


e  =  3c- 2 


With  these  values  it  can  easily  be  verified  by  substitution  into 
equations  (4.22)  that  an  inhomogeneous  steady-state  solution  is: 


Now  we  can  use  conditions  (4.25)  to  test  the  stability  of  this  solution. 
We  obtain  the  following: 


(4.28a) 


(4.28b) 


c  >  1 


Condition  (4.28a)  corresponds  to  (4.25a)  and  (4.28b)  corresponds  to  (4.25b). 
Clearly  by  choosing  c  >  1  we  shall  guarantee  the  stability  of  this  inhomo¬ 
geneous  steady-state  solution. 

Using  the  values  specified  in  equations  (4.26)  let  us  now  look  at 
the  values  of  the  corresponding  homogeneous  steady-state  solution.  We  have: 


X 


0 


2 


(4.29a) 


X 


(4.29b) 


0+ 


2c 


X 


0- 


2c 


(4.29c) 


t 


\ 


■ 

' 
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It  can  easily  be  shown  that  both  X^+  and  are  positive  valued  for 

all  c  >  1  .  Thus  for  c  >  1  ,  all  three  homogeneous  steady-state  solu¬ 
tions  are  physically  meaningful.  By  the  previous  work  that  we  did,  we 
know  that  two  of  them  will  be  stable.  It  turns  out  that  the  solutions 
X^  and  Xq_  are  both  stable  for  all  c  >  1  ,  while  the  X^  =  2  solu¬ 
tion  is  unstable.  Thus  we  see  that  for  this  simple  mechanism  considered 
in  the  two  box  model  it  is  possible  to  have  a  stable  inhomogeneous  steady- 
state  solution  simultaneously  with  a  stable  homogeneous  steady-state 
solution.  Indeed,  because  of  the  existence  of  multiple  homogeneous  steady- 
state  solutions  we  have  at  least  three  simultaneously  stable  solutions  -  two 
homogeneous  and  at  least  one  inhomogeneous  solution. 
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CHAPTER  5 

CONCLUSION 


In  this  paper  we  have  given  evidence  to  support  the  proposi¬ 
tion  that  a  chemical  system  may  under  some  circumstances  exhibit  spatial 
structure.  Recently  there  has  been  reported  an  experimental  confirmation 

g 

of  this  in  the  literature.  Herschkowitz-Kaufman  prepared  a  mixture  of 
Ce^CSO^)^  ,  KBrO^  ,  C^CCOOH^  ,  H^SO^  ,  and  drops  of  Ferroine  as  a  redox 
indicator.  The  mixture  was  made  initially  homogeneous  by  stirring  and 
throughout  the  experiment  the  temperature  was  maintained  constant  at  21 °C  . 
When  the  stirring  was  stopped,  the  solution  immediately  showed  oscillations 
in  time;  the  color  of  the  solution  changed  periodically  from  red  to  blue. 
After  a  number  of  time  oscillations  an  inhomogeneity  in  the  concentration 
appeared  and  gradually  the  solution  arranged  itself  into  alternating  hori¬ 
zontal  layers  of  red  and  blue.  After  about  thirty  minutes  the  system 
approached  chemical  equilibrium  and  again  became  homogeneous. 

Several  things  should  be  noted  about  this  experiment.  An  important 
one  is  that  the  spatial  structure  formed  would  tend  to  support  the  applica¬ 
bility  of  the  derivative  boundary  conditions,  not  the  fixed  boundary  condi¬ 
tions.  For  fixed  boundary  conditions,  the  concentrations  of  the  various 
species  would  have  had  to  remain  constant  at  the  walls  of  the  container. 
Alternating  layers  which  extend  to  the  boundary  preclude  this  possibility. 

On  the  other  hand,  the  horizontal  layers  indicate  that  the  normal  derivative 


of  the  concentration  could  be  zero. 


\ 
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Further  it  is  apparent  from  the  description  of  the  experiment, 
that  the  structure  was  not  caused  by  an  instability  in  the  steady-state 
solution  with  respect  to  spatial  fluctuations.  If  this  were  the  case  the 
structure  would  have  appeared  immediately  after  the  cessation  of  stirring, 
instead  of  being  preceded  by  a  period  of  temporal  oscillations.  It  is 
possible  that  the  steady-state  solution  was  unstable  initially  with  respect 
to  homogeneous  fluctuations  and  so  evolved  into  a  time  dependent  solution 
which  then  became  unstable  with  respect  to  inhomogeneous  fluctuations.  It 
could  also  be  that  this  experiment  is  an  example  of  the  phenomena  discussed 
in  the  fourth  chapter.  That  is,  no  instability  at  all  is  involved,  but 
there  is  simply  a  structured  state  which  is  stable  and  with  the  given  initial 
conditions  the  system  evolved  into  this  state.  Obviously,  further  experimen¬ 
tal  and  theoretical  work  is  needed  to  decide  which,  if  any,  of  these  possibil¬ 
ities  is  correct  in  order  to  explain  the  appearance  of  this  dissipative 
structure. 


In  conclusion,  it  should  be  mentioned  that  it  seems  that  this  area 
of  research  could  be  of  great  interest  in  biology.  It  has  been  shown  that 
the  type  of  mechanism  which  can  lead  to  instabilities  is  consistant  with  the 
mechanisms  established  for  several  biological  processes1'7.  Since  the  biolog¬ 
ical  world  seems  governed  by  structure  (for  example,  the  cell  and  organized 
groups  of  cells),  this  could  provide  a  clue  in  explaining  the  appearance  and 

maintenance  of  this  structure.  Already  the  concept  of  dissipative  instability 

3 

has  been  used  in  a  model  of  a  membrane  .  It  is  to  be  expected  that  further 
applications  will  appear  in  the  literature. 


W 
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